First principle calculations of thermodynamic properties of pure graphene sheet and graphene sheets with Si, Ge, Fe, and Co impurities
Kheyri A1, Nourbakhsh Z2, †,
Plasma Physics Research Center, Science and Research Branch, Islamic Azad University, Tehran, Iran
Physics Department, Faculty of Science, University of Isfahan, Isfahan, Iran

 

† Corresponding author. E-mail: z.nourbakhsh@sci.ui.ac.ir

Abstract
Abstract

The thermal properties of pure graphene and graphene–impurity (impurity = Fe, Co, Si, and Ge) sheets have been investigated at various pressures (0–7 GPa) and temperatures (0–900 K). Some basic thermodynamic quantities such as bulk modulus, coefficient of volume thermal expansion, heat capacities at constant pressure and constant volume of these sheets as a function of temperature and pressure are discussed. Furthermore, the effect of the impurity density and tensile strain on the thermodynamic properties of these sheets are investigated. All of these calculations are performed based on the density functional theory and full quasi harmonic approximation.

1. Introduction

Graphene is a monatomic layer of carbon atoms in a honeycomb lattice. Figure 1 shows the bond length and the σ bonds between the carbon atoms in graphene (C–C bond). The equilibrium bond length between two adjacent atoms in graphene is 0.142 nm.[1] As an element of the graphite structure, graphene has been under investigation since 1916 and the first theory to explain the electronic properties of graphite was formulated in 1947.[2,3] Attempts to exfoliate graphene from graphite began from 1990 while graphene had been considered theoretically for decades. Experimentally, graphene was first separated from graphite in 2004.[4] Graphene has become one of the most famous substances due to its electronic, thermal, optical, and mechanical properties. Graphene is the thinnest and lightest known material, but the characterization of it from the surface science and boundary thermo dynamics looks rather poor. So the purpose of this paper is to investigate the thermal properties of pure graphene sheet and graphene sheets with Si, Ge, Fe, and Co impurities based on density functional theory using the quasi harmonic model.[5]

2. Computational details
2.1. First-principles calculations

The first important conceptual advantage introduced by density functional theory (DFT)[6] is the possibility to describe the ground state properties of a real system in terms of its ground state electronic charge density instead of the far more complicated wave functions. All of the calculations in this paper are performed based on density functional theory using the generalized gradient approximation (GGA) (for the exchange–correlation energy)[7,8] and WIEN2K code.[9] The WIEN2K computational code allows performing electronic structure calculations of solids using density functional theory. It is based on the full-potential linearized augmented plane wave (LAPW) + local orbitals (lo) method, one of the most accurate schemes for band structure calculations. The 4 × 4 × 1 and 2 × 2 × 1 hexagonal honeycomb lattice graphene sheets with 22 Bohr distances between the layers are simulated. By replacing one of the C atoms in the 4 × 4 × 1 and 2 × 2 × 1 graphene sheets with Si, Ge, Fe, and Co impurities, the Si–graphene, Ge–graphene, Fe–graphene, and Co–graphene sheets are created with two different impurity densities (1/32 = 3.1% and 1/8 = 12.5%) (see Fig. 1). In the LAPW + lo approach, each unit cell is divided into two regions: the muffin tin sphere around each atom and the interstitial region outside the spheres.[10,11] The radii of the non overlapping muffin tin spheres are chosen as Rc = 1.13 a.u. and RFe,Co,Si,Ge = 1.5 a.u. The Kmax = 8.5/RMT (a.u.)−1 (RMT is the smallest muffin tin radius and Kmax is the largest augmented plane wave vector) and 180 k-points in the irreducible wedge of the Brillouin zone are chosen. The largest vector of charge density in the potential Fourier expansion in the interstitial region is confined to Gmax = 14.5(Ry)1/2. All of the proposed structures are fully relaxed such that all forces are less than 1 mRy/a.u.

Fig. 1. The (a) 4 × 4 × 1 and (b) 2 × 2 × 1 hexagonal honeycomb lattices of graphene sheets. One atom of carbon in a pure graphene sheet is replaced with Si, Ge, Fe, and Co impurities.
2.2. Quasi harmonic approximation

In the harmonic model, the vibrations of atoms in a crystal are treated as 3nN non-interacting phonons with volume independent frequencies ωi, where n and N are the number of atoms per primitive cell and the number of cells in the solid, respectively. The lack of anharmonicity leads to unphysical behaviors,[12] such as zero thermal expansion, infinite thermal conductivity, etc. The harmonic approximation at any crystal geometry (the quasiharmonic approximation (QHA)) is the simplest way of accounting for the anharmonic effects.[1316] In QHA, the non-equilibrium Helmholtz free energy is

where kB is the Boltzmann constant and ωj is the vibration frequency. The Brich–Murnaghan equation[17] of state is used to fit the calculated energy–volume results

where B0 and are the static bulk modulus and its derivative , respectively, and x = V/V0. Some thermodynamic properties such as constant-volume heat capacity CV, isothermal bulk modulus BT, Grüneisen ratio γth, thermal expansion coefficient α, constant-heat capacity CP, and adiabatic bulk modulus BS are derived directly from F* as

where U = F+TS, S = −(∂F/∂T)V, and .

The full quasi harmonic approximation as implemented in the GIBBS code[5] is used to investigate the thermodynamic properties of pure graphene, Si–graphene, Ge–graphene, Fe–graphene, and Co–graphene sheets with different impurity densities.

3. Results and discussion
3.1. Static properties

To investigate the structural properties of pure graphene sheet and graphene sheets with Si, Ge, Fe, and Co impurities, the total energy as a function of volume for these sheets within GGA is calculated and the results are fitted with the Murnaghan equation of states.[17] The results of these calculations for the impurity density of 3.1% are shown in Fig. 2. The corresponding equilibrium volume per formula unit and bulk modulus (B) of these sheets at T = 0 K are calculated and given in Table 1. The calculated results show that the equilibrium volume of the graphene sheet increases when one of the C atoms in the sheet is replaced by Si, Ge, Fe, and Co impurities. Also the equilibrium volume of these sheets increases with increasing impurity density.

The dominant effect on the bulk moduli of these sheets is from the degree of partial ionicity and strong covalency character. The bulk modulus generally increases with increasing covalency. The ionicity effect is to reduce the bonding charge amount and hence to reduce the bulk modulus.

The bulk modulus of the graphene sheet increases when one of its C atoms is replaced by Co, Fe, and Ge impurities. This is due to the increase of hybridization between the 3d orbitals of the Co, Fe, and Ge atoms and the s, p orbitals of their nearest neighbor C atoms. In our previous paper,[18] it has been shown that the impurity-nearest neighbor distance of the Si–graphene sheet is larger than that of graphene sheets with Co, Fe, and Ge impurities. The bulk modulus of the graphene sheet decreases when one of its C atoms is replaced by a Si impurity. This trend is due to the decrease of covalency with increasing nearest neighbor distance and larger Si atomic radius than C atomic radius.

Furthermore, the bulk moduli of these sheets decrease and are in the same order of magnitude when the impurity density increases to 12.5%. We also apply a small strain along the diagonal of the sheet to investigate the thermodynamic quantities at tensile strain. The mechanical properties and shear modulus of graphene have been reported,[1930] however to the best of our knowledge, there are no experimental or theoretical values of the equilibrium unit cell volume and bulk modulus of this sheet with Si, Ge, Fe, and Co impurities reported.

Table 1.

The calculated equilibrium volume and bulk modulus of pure graphene and impurity–graphene sheets (impurity = Si, Ge, Fe, and Co with densities of 3.1% and 12.5%) at 0 K.

.
Fig. 2. Total energy versus volume curves of (a) pure graphene sheet and (b)–(e) graphene sheets with 3.1% Si, Ge, Fe, and Co impurities.
3.2. Thermal properties

The thermal properties of pure graphene sheet and graphene sheets with Si, Ge, Fe, and Co impurities are investigated in the temperature range of 0–800 K (below the melting point of pure graphene sheet) and the pressure range of 0–7 GPa. The pressure effect on the volume (EOS curves) at room temperature (299 K), 500 K, and 800 K are calculated and the results are shown in Fig. 3. The volume of these sheets decreases rapidly with increasing pressure. At a specific pressure, the volume of these sheets increases with increasing temperature.

A better presentation of these changes is shown in Fig. 4. The volume of these sheets increases with increasing temperature and decreasing pressure but the rate of increase is moderate, and the volume as a function of temperature for all of these sheets is nearly linear at various pressures. The bulk moduli of pure graphene sheet and graphene sheets with Si, Ge, Fe, Co impurities as a function of temperature at various pressures are calculated and the results are shown in Fig. 5. It shows that the bulk moduli of these sheets decrease with increasing temperature at a given pressure with a relatively mild slope and increase with increasing pressure at a given temperature. As mentioned, the bulk modulus is a gauge of the crystal hardness. The electron charge distributions of pure graphene sheet, Si–graphene, Ge–graphene, Fe–graphene, and Co–graphene sheets in the (001) plane at two different pressures (0 GPa and 7 GPa) are demonstrated in Fig. 6. At zero pressure, the electron charge distribution of pure graphene sheet is spherical with bonding around all atoms. When one of the C atoms in the graphene sheet is replaced by Si, Ge, Fe, and Co impurities, the spherical symmetry of the electron charge distribution around the impurities decreases due to the electro negativity difference between C (2.55) and impurity atoms (Fe (1.83), Co (1.88), Si (1.90), and Ge (2.01)). There is a partial ionic and strong covalent character along the impurities and their nearest neighbor C atoms. The charge transfer from the C atom to the impurities increases with increasing pressure, so the bulk modulus of these sheets increases with increasing pressure.

Fig. 3. The relationship between volume and pressure of (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets at 299 K, 500 K, and 800 K.

The heat capacities of pure graphene, Si–graphene, Ge–graphene, Fe–graphene, and Co–graphene sheets at constant volume (CV) and constant pressure (CP) as a function of pressure are investigated and shown in Fig. 7. As shown in Fig. 7, the behavior of CV is exactly like that of CP for all these sheets. CV and CP increase with increasing temperature and decrease slowly with increasing pressure and drop rapidly at room temperature. At constant pressure and temperature, CV and CP of these sheets increase when the impurity density increases from 3.1% to 12.5%. As can be seen, increasing the impurity density has no significant effect on CP at high temperatures.

Fig. 4. The relationship between volume and temperature of (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets at different pressures.
Fig. 5. The relationship between bulk modulus and temperature of (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets at 0 GPa, 3 GPa, 5 GPa, and 7 GPa.
Fig. 6. The electron charge density of (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets in the [001] plane for two different pressures of 7 GPa (right) and 0 GPa (left).

The CV and CP of pure graphene, Si–graphene, Ge–graphene, Fe–graphene, and Co–graphene sheets at tensile strain from 0 to 4 GPa are calculated and compared with the corresponding results under hydrostatic pressure in Fig. 8. The results show that CV and CP of these sheets under tensile strain and hydrostatic pressure have similar behaviors. The CV and CP of these sheets under hydrostatic pressure are larger than the corresponding values of these sheets under tensile strain at low temperature. At definite pressure, the CV and CP of these sheets under tensile strain tend to the CV and CP of these sheets under hydrostatic pressure with increasing temperature.

The CV of these sheets with different impurity densities as a function of temperature at different pressures is plotted in Fig. 9. As expected, for all these sheets, CV is proportional to T3 at low temperature and tends to the Dulong and Petit limit at high temperature (around 780 J/mol·K). The Dulong and Petit limit increases with increasing impurity density.

Fig. 7. The heat capacities at constant volume (CV) and constant pressure (CP) as a function of pressure at different temperatures for (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene and (e) Co–graphene sheets. The black lines are for 3.1% impurity density and the blue ones are for 12.5%.
Fig. 8. The heat capacities at constant volume (CV) and constant pressure (CP) as a function of pressure at different temperatures for (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets.
Fig. 9. The heat capacity at constant volume (CV) as a function of temperature for (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets at different pressures. The black lines are for 3.1% impurity density and the blue ones for 12.5%.
Fig. 10. The coefficient of volume thermal expansion (α) as a function of temperature for (a) pure graphene, (b) Si–graphene, (c) Ge–graphene, (d) Fe–graphene, and (e) Co–graphene sheets at different pressures.

The Dulong and Petit limit of these sheets decreases by about 1.18% with increasing pressure. The CV of these sheets increases with increasing impurity density at given pressure and temperature. The CV of these sheets decreases slightly with increasing pressure. The coefficient of volume thermal expansion (α) for pure graphene and graphene sheets with Si, Ge, Fe, Co impurities as a function of temperature at different pressure is calculated and shown in Fig. 10. The volume thermal expansion coefficient arises from the anharmonicity of the interatomic forces. As can be seen from Fig. 10, for a given pressure, at low temperature (between 100 K and 200 K), α increases rapidly, and at high temperature, it tends to increase linearly. Also α becomes smaller with increasing pressure for a given temperature.

4. Conclusion

The volume of pure graphene sheet and graphene sheets with Si, Ge, Fe, and Co impurities increases with increasing temperature at a given pressure. The bulk modulus of pure graphene sheet at T = 0 K increases when one of the C atoms in the graphene sheet is replaced by Ge, Fe, and Co impurities and decreases when it is replaced by a Si impurity. The bulk modulus of these sheets decreases with impurity density increasing. The bulk modulus of these sheets under tensile strain increases compare to the corresponding value under hydrostatic pressure.

With increasing pressure and temperature, the bulk modulus of these sheets decreases and tends to a constant value. CV is proportional to T3 at low temperature and tends to the Dulong and Petit limit at high temperature. CV and CP of the strained sheets significantly decrease compare to those of these sheets under hydrostatic pressure at low temperature. α increases rapidly at low temperature (between 100 K and about 200 K) and tends to increase linearly at high temperatures for a given pressure.

Reference
1Heyrovska R2008arXiv:0804.4086
2Novoselov S2010Scientific Background on the Nobel Prize in Physics 2010, GRAPHENEThe Royal Swedish Academy of Science
3Wallace P R 1947 Phys. Rev. 71 622
4Novoselov K SGeim A KMorozov S VJiang DZhang YDubonos S VGrigorieva I VFirsov A A 2004 Science 306 666
5Otero-de-la-Roza AAbbasi-Pérez DLuaña V 2011 Comput. Phys. Commun. 182 2232
6Kohn WSham L J 1965 Phys. Rev. 140 A1133
7Perdew J PBurke KErnzerhof M 1996 Phys. Rev. Lett. 77 3865
8Perdew J PBurke KErnzerhof M1997Phys. Rev. Lett.781396
9Blaha PSchwarz KMadsen G K HKvasnicka DLuitz J2014WIEN2k, An Augmented Plane Wave plus Local Orbitals Program for Calculating Crystal PropertiesViennaUniversity of Technology
10Anderson O K 1975 Phys. Rev. 12 3060
11Loucks T LBenjamin W A1967Augmented Plan Wave MethodNew YorkW.A. Benjamin, Inc
12Ashcroft N WMermin N D1976Solid State PhysicsNew YorkThomson Learning Inc
13Chen X JZhang CMeng YZhang R QLin H QStruzhkin V VMao H K 2011 Phys. Rev. Lett. 106 135502
14Umemoto KWentzcovitch R MSaito SMiyake T 2010 Phys. Rev. Lett. 104 125504
15Wood B CMarzari N 2009 Phys. Rev. Lett. 103 185901
16Oganov A RChen JGatti CMa YGlass C WLiu ZYu TKurakevych O OSolozhenko V L 2009 Nature 457 863
17Birch F 1947 Phys. Rev. 71 809
18Keyri ANourbakhsh ZDarabi E 2016 J. Supercond. Nov. Magn. 29 985
19Lee CWei XKysar J WHone J 2008 Science 321 385
20Frank I WTanenbaum D MVan der Zande A MMcEuen P L 2007 Vac. Sci. Technol. 25 2558
21Liu FMing PLi J 2007 Phys. Rev. 76 064120
22Lier G VAlsenoy C VDoren V VGeerlings P 2000 Chem. Phys. Lett. 326 181
23Zhao HMin KAluru N R 2009 Nano Lett. 9 3012
24Zhao HAluru N R 2010 J. Appl. Phys. 108 064321
25Sakhaee-Pour A 2009 Solid State Commun. 149 91
26Zakharchenko K VKatsnelson M IFasolino A 2009 Phys. Rev. Lett. 102 046808
27Tsai JTu J 2010 Mater. Des. 31 194
28McSkimin H JAndreatch P 1972 J. Appl. Phys. 43 2944
29Roundy DCohen M L 2001 Phys. Rev. 64 212103
30Zhang YSun HChen C 2004 Phys. Rev. Lett. 93 195504